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Abstract 


The stellar initial mass function (IMF) impacts nearly all observable 
properties of galaxies [1], controls the production rate of heavy ele- 
ments, and governs how much energy is available to regulate galaxy 
growth. Theoretical work predicts that the high-redshift IMF may be 
more top-heavy compared to the local Universe, due to higher gas 
pressures, a higher Cosmic Microwave Background temperature, and 
lower metallicities [2, 3]. However, direct observational evidence for a 
top-heavy IMF at high-redshift remains elusive. Here we report the 
detection of two Lyman-a-emitting galaxies at redshift 5.9 and 7.9 
that show evidence for exceptionally top-heavy IMFs. Our analysis of 
JWST/NIRSpec data demonstrates that these galaxies exhibit spectra 
which are completely dominated by the nebular continuum. Alongside 
a clear Balmer jump, we observe a steep turnover in the ultraviolet 
continuum. Although this feature can be produced by an extremely 
thick damped Lyman-a system with holes, we show instead that this 
turnover is two-photon emission from neutral hydrogen. Two-photon 
emission can only dominate if the ionizing emissivity is > 10x that ofa 


typical star-forming galaxy. While weak He II emission disfavours ion- 
izing contributions from AGN or X-ray binaries, such radiation fields 
can be produced in star clusters dominated by low-metallicity stars of 
2 50 Mo, where the IMF is 10 — 30x more top-heavy than typically 
assumed. Such a top-heavy IMF implies our understanding of star for- 
mation in the early Universe and the sources of reionization may need 
revision. 


The detection of numerous Lya-emitting galaxies at z > 5.5 from the JWST Advanced 
Deep Extragalactic Survey (JADES) was recently reported in [4]. Among this sample 
we identify one object, JADES-GS+53.12175-27.79763 (GS-NDG-9422 hereafter), at 
z = 5.943, that exhibits a clear spectral discontinuity (Figure 1) of 15.0 + 0.9 nJy 
(observed frame) at the location of the Balmer jump (Arest = 3645 A). This feature 
arises due to electron recombinations with ionized hydrogen to the first excited state 
and only appears in spectra of galaxies that contain young stellar populations with 
high ionizing photon production efficiencies (ion). Balmer jumps have been detected 
in galaxies at low-redshift [5, 6] and are predicted by numerical simulations [7] and 
spectral energy distribution (SED) fitting of JWST photometry [8] to be common 
at high-redshift. In contrast, the other nebular continuum components, free-free and 
two-photon emission, are predicted to be subdominant compared to the stellar and 
free-bound contributions in high-redshift spectra. 

The subdominance of two-photon emission arises because the Zon values of typi- 
cal low-metallicity stellar populations are not high enough to overcome the steep UV 
stellar continuum slopes [9]. In contrast, very hot stars, with blackbody temperatures 
of ~ 100,000 K, produce enough ionizing photons to power nebular continuum emis- 
sion that outshines the stellar UV emission, and are predicted to exhibit a strong 
two-photon continuum bump, peaking at Arest © 1430 A [10-12]. 

In addition to the strong Balmer jump, GS-NDG-9422 exhibits a steep turnover in 
the continuum at Arest zi 1430 A (Figure 2). Similar UV turnovers are often observed 
as a result of absorption from foreground neutral hydrogen — either the neutral 
intergalactic medium (IGM) [13] or a Damped Lyman-a (DLA) system [14]. Fitting 
GS-NDG-9422 with a typical stellar population requires excess attenuation beyond 
IGM damping (Figure 2) and necessitates DLA column densities > 107° cm~?, higher 
than any known DLA. Furthermore, such high column densities imply zero trans- 
mission at 1216 A, conflicting with the strong observed Lya emission. This could be 
reconciled by invoking a DLA with 30 % leakage (Methods), however the plausibility 
of such an extreme column density with a low covering fraction is unclear. Further- 
more, explaining this with a fine-tuned geometry would seemingly be in tension with 
the identification of two other objects, the Lynx arc at z = 3.36 [15] and a recently- 
observed galaxy (A2744-NDG-ZD4) from the UNCOVER program [16] at z = 7.88, 
that exhibit a similar turnover in the UV continuum and strong Lya (Figure 2). 

Rather, we argue that, like the Balmer Jump, the observed UV turnover is of 
nebular origin and a signature of the two-photon continuum arising from 2s > 1s 
transitions in neutral hydrogen. To demonstrate this, the top panel of Figure 3 shows 
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Fig. 1 2D (top) and 1D (middle) spectra of GS-NDG-9422 showing clear evidence of a spectral 
break at 2.53 wm. The grey shaded region shows noise in the observed spectrum, which is well below 
the signal at the location of the break. Bottom: model nebular continuum at 18,300 K showing the 
relative contribution of the two-photon, free-bound, and free-free components. This model continuum 
is also shown as the purple dashed line in the middle panel, renormalized to the observed spectrum. 


the predicted nebular continuum for gas at a density of 100 cm~? at a range of temper- 
atures. We highlight in magenta the prediction for gas at a temperature of 18, 300 K, 
equal to that measured from [O m1] A4363/A5007. The fact that the predicted con- 
tinuum from an independent temperature measurement nearly perfectly matches the 
observed continuum indicates that the UV turnover is undoubtedly a detection of the 
two-photon nebular continuum: the spectrum of GS-NDG-9422 is entirely dominated 
by nebular emission. 

In order for both the two-photon and free-bound continua to dominate over the 
stellar spectrum in the rest-frame UV and optical, the source population must be 
emitting significantly more ionizing photons than standard SSP models, necessitating 
blackbody temperatures of T = 90,000 K (see Methods), much hotter than massive O 
stars (up to ~50,000 K; [17]). Hence, if the ionizing sources illuminating these three 
high-redshift galaxies are of stellar origin, they must be exotic in nature. 
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Fig. 2 Top. Zoom in on the rest-frame UV region of the JWST NIRSpec PRISM spectrum of GS- 
NDG-9422 (black) compared with the best fit models using a standard SSP (magenta) accounting 
for both the z = 6 IGM opacity (dashed yellow) and a DLA with column density of 1029! cm7? 
(dotted cyan). The full spectral fit for these models is shown in Methods. Bottom. Zoom in on the 
rest-frame UV region for GS-NDG-9422 (black) compared with the spectra of the z = 3.35 Lynx arc 
(left) and A2744-NDG-ZD4 at z = 7.88 (right). We have normalised the spectra and shifted them all 
to the rest frame. 


Spectral fitting of GS-NDG-9422 with photoionization models (second panel of 
Figure 3) demonstrates that a blackbody model with T = 10505 K can provide a 
good fit to both the shape of the continuum and the majority of lower-ionization 
state emission line ratios. However, where our simple blackbody model fails is that it 
strongly overpredicts the flux of He 11 lines compared to those observed in the spectrum 
(bottom panels of Figure 3). The weak He 11 A1640 and He 11 \4686 in GS-NDG-9422 
places tight constraints on the ionizing spectrum at E > 4 Rydberg. To explore this, 
we run models varying the leakage fraction of photons with E > 4 Rydberg from 0% 
(i.e. no He*-ionizing photons) to 100% (unattenuated blackbody). We conclude that 
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Fig. 3 First Panel. Predicted nebular continuum (coloured lines) for gas of primordial composition 
at fixed temperature (as listed in the colour bar and legend) compared to the JWST NIRSpec PRISM 
spectrum of GS-NDG-9422. The magenta line shows the prediction when using the gas temperature 
(To3) derived from the ratio of [O m1] 44363/\5007, which nearly perfectly matches the continuum 
of the observed spectrum. Second Panel. Spectrum of GS-NDG-9422 (black) compared with the 
best fit blackbody model (magenta). The blue/green lines show the blackbody models with different 
leakage fractions of He 11 ionizing flux (see Methods). Third Panel. Spectrum of GS-NDG-9422 (black) 
compared with the Wolf-Rayet model with the most similar spectrum to the inferred SED (dashed 
orange) or various individual massive Pop. III stars with different effective temperatures (dotted 
lines). Fourth Panel. Zoom in on the regions around He 11 \1640 and He 11 \1486 demonstrating why 


the He 1 ionizing flux needs to be reduced in the blackbody models. Lines are the same as in the 
second and third panels. 
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Fig. 4 Comparison of the best fit blackbody SED model (black) with the predicted spectra of 
Pop. III stars with varying IMFs (top left), individual massive Pop. III stars with different effective 
temperatures as shown in the colour bar (top right), low-metallicity Wolf-Rayet stars (bottom left), 
and stripped stars (bottom right). The dotted black line shows the blackbody spectrum prior to 
extinguishing the He 1 ionizing radiation. 


75%—80% of the Het-ionizing photons emitted by the blackbody must be extinguished 
to match the observed He 1 emission (i.e. leakage fraction of S25%; Figure 4). 

We note that other sources, including active galactic nuclei (AGN) or high-mass 
X-ray binaries (HMXBs), can have significant high-energy photon outputs, and have 
been invoked to explain peculiar emission line ratios seen at high-redshift [18, 19]. 
However, the weak He I emission disfavours power-law SEDs that extend past the 
He II ionizing edge. We exclude HMXBs due to the strong over-prediction of He 1 
emission in these models (see Methods). Reproducing the observed spectrum with an 
AGN would require a peculiar absorption of the power-law spectrum [20]. Furthermore, 
the emission lines in GS-NDG-9422 do not show any evidence of a broad component 
that would be indicative of a broad-line region. 

Massive metal-free Population III stars are predicted to have sufficiently high 
temperatures to power a two-photon dominated spectrum [10]. While the IMF of Pop- 
ulation III stars is highly uncertain [21], comparing extremely top-heavy (Pop. III.1) 
and moderately top-heavy (Pop. II.2) Yggdrasil models [12] with our extinguished 
blackbody (Figure 4 top left) indicates these produce too many Het ionizing pho- 
tons, while predictions for even less top-heavy IMFs begin to fall short of the required 
ĉion: In contrast, some individual Pop. III star models from [22] with effective tem- 
peratures of ~ 97,000 K and masses of 85 Mọ to 108 Mo reasonably reproduce the 
required ionizing SED (Figure 4 top right). We note that the measured metallicity of 
12+ log,)(O/H) = 7.59+ 0.01 suggests the presence of Pop. III stars is unlikely, unless 
we are witnessing the immediate enrichment and illumination of metals produced by 
primordial stars. Nevertheless, we consider these Pop. III star models in our analysis 


under the assumption that the atmospheres may be generally representative of hot 
massive stars at low metallicity. 

Wolf-Rayet stars or stripped stars not only exhibit surface temperatures above 
10° K, but can also have helium atmospheres that provide the necessary opacity to 
reduce their He" ionizing output. Compared to our best-fit extinguished blackbody, 
we find that stripped star SEDs [23] produce too few He" ionizing photons. How- 
ever, select low-metallicity Wolf-Rayet models [24] show a remarkable resemblance 
(Figure 4). 

We note that spectral fitting of some low-redshift reionization analogues has sug- 
gested a similar need for significant ionizing contribution from a hot (T > 80,000 K) 
blackbody [25], however the key difference in the high-redshift galaxies presented here 
is the dominance of the two-photon continuum. This necessitates an extremely high 
ĉion that cannot be produced if a substantial fraction of the ionizing photons are 
emitted by typical OB stars. 

The progenitor masses of Wolf-Rayet stars at the metallicity of GS-NDG-9422 are 
not well constrained [26]. Masses of > 37 Mo have been estimated for Wolf-Rayet stars 
in the Small Magellanic Cloud (SMC) with T = 100, 000 K [27], implying even higher 
progenitor masses. While this search may not be exhaustive, we conclude from Figure 4 
that the nebular-dominated spectrum of GS-NDG-9422 is consistent with ionization 
powered by low-metallicity massive stars (= 50 Ma), perhaps in the Wolf-Rayet phase. 

Assuming a typical IMF, we expect to form one such 50 Mo star per ~1,300 Mo of 
stellar mass. We repeat our photoionization modelling with our best-fit low-metallicity 
massive star and Wolf-Rayet models, maximizing the mass of an underlying stellar 
population with a normal IMF such that model spectrum remains a good fit (see 
Methods). For the two-photon continuum to remain dominant, we require < 112 Mo 
of a normal stellar population per one low-metallicity massive star, or < 137 Mo 
per one Wolf-Rayet star. This necessitates 35 times, or 9.5 times, more massive stars 
compared to a typical IMF in the low-metallicity massive star scenario, or Wolf-Rayet 
star scenario, respectively (Methods). Either scenario implies that the IMF in GS- 
NDG-9422 must be extraordinarily top-heavy. 

The IMF is predicted to get progressively more top-heavy for low-metallicity gas 
at high pressure [2, 3], and a top-heavy IMF has been suggested to explain IR emission 
line ratios [28] and the surprising abundance of UV-bright galaxies [29] observed at 
high redshift. 

The strong contribution of the two-photon continuum to the overall spectrum of 
GS-NDG-9422 demonstrates the need to update SED fitting templates to account 
for stellar populations with extreme Zon, The intrinsic Zon for our best-fit black- 
body model is 26.4 Hz erg), 10x greater than commonly adopted for star-forming 
galaxies in reionization models [30]. Since our modelling indicates that GS-NDG- 
9422 has an ionizing escape fraction of 7.3%, significantly fewer nebular dominated 
galaxies (NDGs) are required to drive reionization compared to star-forming galaxies, 
motivating future work to constrain their number density and luminosity function. 

The striking similarities between GS-NDG-9422 and both the Lynx arc and A2744- 
NDG-ZD4 suggest that NDGs powered by a very top-heavy IMF may not be rare 
at high redshift. Selecting NDGs at high redshift when the IGM is opaque to Lya 


will be challenging because without Lya, the two-photon UV turnover is degenerate 
with a DLA. Nevertheless, the exotic nature of GS-NDG-9422 and A2744-NDG-ZD4 
demonstrate that our understanding of star formation in the early Universe requires 
substantial revision. 


1 Methods 


1.1 Data 


GS-NDG-9422 was observed in the JADES survey (PID: 1210, PI: Luetzendorf) in five 
spectral modes, receiving 28 hours integration in Prism/CLEAR and 7 hours integra- 
tion in each of G140M/FO70LP, G235/F170LP, G395M/F290LP and G395H/F290LP. 
We use the reduced spectra released as part of the JADES Public Data Release [31, 32]. 

A2744-NDG-ZD4 was observed as part of the Ultra-deep NIRCam and NIRSpec 
Observations Before the Epoch of Reionization (UNCOVER) program (ID: 2561, PI: 
Labbe; [33]) in the lensing field Abell 2744. A2744-NDG-ZD4 is lensed with magni- 
fication factor of u = 1.96 [34], and received 2.3 hours inegration in Prism/CLEAR. 
Reduced data was retrieved from the Dawn JWST Archive (DJA). The data was 
reduced using the custom-made pipeline MsaExp v.0.6.12 7. Further details of the 
reduction pipeline are given in [35]. 

The spectrum of the Lynx arc was taken with the Keck I telescope using the Low 
Resolution Imaging Spectrograph (LRIS) [15]. 


1.2 Emission line fitting 


Where possible, emission lines were fit with a single component Gaussian profile with 
the continuum modelled as a 1D spline. In cases where lines are sufficiently blended 
we fit the whole complex with one component. In some cases, partially blended lines 
are fit simultaneously with neighbouring lines and fluxes reported separately. These 
are marked in Table 2. We fit all identifiable lines and report upper limits for notable 
undetected lines. 

Line fluxes from higher resolution grating spectra of GS-NDG-9422 were measured 
independently. Fluxes derived from different observations are mildly discrepant in some 
cases. Notably, H6, [Om] A4959, [Orm] 45007 and Ha lines exhibit higher fluxes in 
the grating modes. This behaviour is reported in [32] who suggest that the Prism flux 
calibration is more reliable. We note that in low-resolution data, the continuum level 
for some emission lines can be ambiguous, especially for Lya, He 1 + O m1] and [O 1], 
which introduces uncertainty into the emission lines flux. Throughout our analysis, we 
adopt the Prism fluxes where possible. However, we ensure that conclusions presented 
in this work are generally consistent with the measured grating fluxes. Emission lines 
in the high-resolution grating show no evidence of a broad component and are well fit 
with a single component with velocity dispersions < 200 km s+. 


https: //dawn-cph.github.io/dja/index.html 
*https: //zenodo.org/record/8319596 


Prism/CLEAR G140M G235M G395M G395H 


Ly-a 110.1 £2.8 101.4 4.3 
Nv \1239 < 5.4 
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Niv] A1483 < 2.4 
Niv] A1486 < 3.1 

Civ 41548 21.5+1.17 

Crv 41550 17.6 + 1.17 
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Table 1 Emission line flux measurements for GS-NDG-9422 across each observed 


spectrum. Fluxes are reported in units of x10~!9 erg s~! cm~?. 


+ Fit simultaneously with neighbouring line and fluxes reported separately. 
t Blended with He 1 A4713. 


1.3 Calculation of physical conditions 


Throughout this section, we make use of PYNEB [36] using atomic data from CHIANTI 
(version 10.0.2; [37, 38]). 


Dust extinction 


Balmer decrements Hô/H8 and Ha /H from the Prism and H9/H{ from the grating 
are consistent with Case B values, indicating that there is no significant dust reddening 
in GS-NDG-9422 (Figure 5). Measured Ha/Hf ratios are lower than theoretically 
predicted at T = 104 K. Note, this is not suggestive of dust reddening, since this would 
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Fig. 5 Observed Balmer decrement ratios compared to theoretical values at different temperatures 
and densities. Green squares give the measured ratios from the Prism, while diamonds give ratios 
measured from gratings, where available. Solid black lines give the theoretical ratio for case B recom- 
bination at Te = 104 K and ne = 100 cm™3. Coloured points show theoretical ratios for a range of 
temperatures (104 < Te < 2.5 x 104 K; color) and densities (ne = 10, 100, 1000 cm~?; marker size). 


act in the opposite direction. At higher temperatures, the theoretical ratio decreases, 
and our G395H/F290LP measurement is consistent with theoretical ratios with Te Z 
2x 10t K, possibly indicative of a very hot nebula. Hy/H@ measured from the medium- 
resolution grating is marginally below the theoretical value. In isolation, this could 
suggest non-zero dust reddening; however, this evidence is outweighed by the other 
measured ratios. The high-resolution grating returns a much lower Hy/H@ = 0.3+0.03. 
However, we note that the continuum is undetected in the high-resolution grating, 
which contributes uncertainty to the measured ratio. Note that much of the focus of 
this work is the unusually strong two-photon nebular continuum at Arest ~ 1400 Å. If 
there is significant dust in this system, the intrinsic two-photon continuum would be 
even higher than measured here, requiring more extreme scenarios than those discussed 
in this work. 


Electron temperature 


The temperature-sensitive [O 111] 4363/5007 ratio can be measured from each of the 
Prism, G395M and G395H observations, yielding three consistent independent temper- 
ature measurements (Te = 1.83-40.15, 1.99+0.18, and 1.81-40.18x 10t K, respectively). 
The temperature derived from the medium-resolution [O 1m1] A1666/A5007 ratio is 
somewhat lower (Te = 1.70*(¢ x 104 K). However, the He 1+O 1m] flux mea- 
sured from the medium-resolution grating is significantly lower than that of the 
Prism. Instead, the measured temperature from above (Te = 1.83 x 104 K) implies 
AX1660,1666/A5007 = 0.08, which gives f\,1660,1666 = 8.2 + 0.1 x 10719 erg s~+ om? 
based on the [O 1m1] 45007 Prism flux, consistent with O 111] contributing ~55 % of 
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Property Derived value 


Te(OTT) 18, 300 £ 1,500 K 
Ne SS 10° em? 
12+ log(OTtt+/Ht) 7.58 £0.01 
12+log(O+/H+) 5.794963 

12 + log(O/H) 7.59 + 0.01 
log(NT* /OTF) < —0.85 
log(Ct+ /O+4 —0.76 + 0.03 
log(C/O) —0.73 + 0.03 
log(Ne/O) —0.8510-03 
Het /Ht 0.10 + 0.01 
Het+ /Ht 0.005 + 0.001 
He/H 0.11 + 0.01 
log(Ot+/Ot+) 1.7970 oa 

Det" /Het 0.06 + 0.01 


Table 2 Physical properties and chemical 
abundances derived for GS-NDG-9422. 


the measured Prism He 11+O 111| blend. The implied He 11 A1640 flux (6.9 + 1.1 x 
10-19 erg al cm7?) is consistent with that expected by extrapolating from fy4686 
using the theoretical He 11 1640/4686 ratio. Note that this further supports the 
conclusion of a lack of dust in this system since the He 11 A1640 would be sub- 
ject to extremely high extinction, relative to the A4686 line. The magnitude of the 
Balmer jump is sensitive to Te (e.g. [6, 39]). In Figure 3, we show this behaviour 
for model nebular continuum spectra for a range of temperatures, finding that the 
Balmer jump magnitude and two-photon continuum contribution are fully consistent 
with Te = 1.83 + 0.15 x 10* K, as derived from [O 1]. 


Electron density 


The detection of the two-photon nebular continuum implies ne < 10? om "7 since this 


feature is strongly suppressed by /-changing collisions at higher densities (Figure 8). 
The density-sensitive doublets of C 111] and [O 11] are not resolved in our observations, 
while N Iv] and [S 11] are not detected. We report a marginal detection of the [Ar rv] 
AA4711, 4740 doublet in the Prism spectrum. These lines are partially blended with 
each other and He 11 4686, while [Ar rv] 44711 is also completely blended with He 1 
A4713. Our three-component fit to this complex yields [Ar Iv] \4711/A4740 = 1.6£0.8 
after subtracting the predicted He 1 4713 contribution (assuming \4713/A4471 = 
0.15), consistent with the low-density limit (ne < 10+ cm~’). 


~N 


Chemical abundances 


We derive chemical abundances for GS-NDG-9422 adopting Te = 1.83 x 104 K and 
a density of ne = 100 cm~® following the procedure in [40]. Given the apparent 
agreement between T.(H*) and T.(OtT), we assume a constant temperature for all 
ionisation zones, deriving 12+log(O/H) = 7.59+0.01 with both Prism and grating line 
fluxes. We derive the carbon abundance from the C 11] AA1907, 1909 / [O ml 45007 
ratio measured from the Prism, assuming no dust, finding log(C/O) = —0.73 + 0.03 
after applying the ionisation correction factor presented in [41]. However, we note the 
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significant emission from C Iv in the spectrum. Given the high Zon implied for GS- 
NDG-9422, the ICF may not be representative of the conditions in GS-NDG-9422 (e.g. 
[42]), and the quoted C/O may represent a lower limit. 

Since no nitrogen lines are detected in GS-NDG-9422, we can only place upper 
limits on the nitrogen abundances. The low OF abundance implies that the N* abun- 
dance is likely negligible. We instead consider the N** abundance using N 111] AA1750 
/ [O 1m] A5007 measured from the Prism, and N 1] AA1750 / O 11] 41666 measured 
from the grating. These yield 3-o upper limits of log(N**/O**) < —0.85 and < —1.01 
respectively. We adopt the former as our preferred limit. 

We detect several helium recombination lines from both the singly and doubly 
ionised states. Prism measurements of He 1 \4471 / HO and He 1 \5875 / HO yield con- 
sistent measurements of Het /H+ = 0.10+0.005 and 0.10+0.01 respectively. Deriving 
a Het*/H*+ abundance using the He 11 4686 line results in only a 6+ 1 % contribu- 
tion to the total helium abundance, giving Hei" /Ht+ = 0.005 + 0.001. Together, this 
implies a total helium abundance of He/H = 0.11 + 0.01, higher than typical values 
observed in low-metallicity systems [43], but consistent with some massive globular 
clusters [44]. However, we caution that our derived helium abundance does not account 
for collisional emission or self-absorption, which could be significant (e.g. [45]). 


1.4 DLA with a Contrived Geometry 


We consider what kind of DLA geometry can reproduce the spectral shape of GS-NDG- 
9422 and A2744-NDG-ZD4. For GS-NDG-9422, we use the best fit BPASS model (see 
below) and apply damping due to a foreground DLA with leakage through optically- 
thin channels. A similar approach is applied to A2744-NDG-ZD4; however, due to 
the lower signal-to-noise, we apply a power-law fit to the spectrum to represent the 
intrinsic stellar continuum. In the top panel of Figure 6 we show the fits to both 
galaxies, and demonstrate that they require column densities > 102°? cm~? and a 
covering fraction of 70%. Not only would these DLAs be the highest neutral columns 
observed (see bottom panel of Figure 6), they must also have optically thin channels 
covering 30% of the total area. Furthermore, at these column densities, the gas can 
self-shield from the local-radiation field, and the core would be expected to be fully 
molecular. Given the fine-tuned requirements, we consider this scenario of a thick DLA 
with holes highly unlikely. 


1.5 Spectral Modelling 
Spectral modelling is performed using the photoionization code CLOUDY v23 [49]. 


Nebular continuum as a function of temperature 

To calculate the nebular continuum in the top panel of Figure 3, we set up a spher- 
ical shell of gas with primordial composition in CLOUDY and irradiate it with two 
“lasers” with energies just above that needed to ionize hydrogen and singly ionize 
helium. The gas temperature is varied in the range 10415 < T/K < 10435, 
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Fig. 6 Top. Fits to GS-NDG-9422 and A2744-NDG-ZD4 using a DLA model with a 70% covering 
fraction. Bottom. Comparison of the required DLA column density with low-redshift star-bursting 
galaxies (orange, [46]), intermediate redshift GRBs (black, [47]), and high-redshift JWST galaxies 
(blue, green, [35, 48]), as well as the galaxies presented in this work (magenta), that also require a 
< 100% DLA covering fraction. 


Models with standard SSPs 


We adopt the BPASS v2.2.1 SSP models including binary stellar evolution [50] with 
a maximum mass of 300 Mo and a high-mass slope of —2.35. We consider a population 
with an age of 3 Myr and a metallicity of 0.1Zọ, consistent with that measured for the 
gas from the spectrum. We adopt a spherical geometry with an inner radius of 0.1 pc. 
The calculation is stopped at an electron fraction of 1% or when the neutral column 
density reaches 10187 cm~?, consistent with the minimal observed Lya emission offset 
(< 100 km s~?) [51]. We vary gas density, ionization parameter, metallicity (assuming 
solar abundance patterns from [52]), and carbon abundance. The gas temperature is 
fixed to that measured from [O 111] 4363/5007. The best fit model (n = 10° em", 


13 


T T 
— GS-NDG-9422 ---- BPASS + IGM 
Best Fit BPASS = BPASS + DLA (1023-1 cm~?) 


10719 


fy ferg/s/cm?/A] 


1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 
Observed wavelength [um] 


Fig. 7 JWST NIRSpec PRISM spectrum of GS-NDG-9422 (black) compared with the best fit models 
using a standard SSP (magenta) accounting for both the z = 6 IGM opacity (dashed yellow) and a 
DLA with column density of 10?3-1 cm~? (dotted cyan). 


logip(U) = 1.2, logy9(Zo/Zo) = —1.1, and logio(Zc/Zo) = —1.4; Figure 2 & 7) is 
obtained matching the line strengths of [O 111] 45007, [O 11] AA3727, C 111] AA1909, Ha, 
and Hy with respect to HS, meaning the resulting continuum is a prediction of the 
model. Assuming z = 6 IGM transmission curves from [53], we find that IGM damping 
is unable to produce the observed UV turnover from this BPASS model (yellow dashed 
line in Figure 2). Even the maximial IGM damping wing, calculated following the 
formalism presented in [13, 54], cannot reproduce this feature. Similar results were 
found in [14] at much higher redshift. We model damping due to the presence of a 
DLA with CLOUDY by calculating the transmission of the best-fit model through 
slabs of neutral hydrogen with varying column densities (best fit shown as dotted blue 
line in Figure 2 & 7). 


Models with blackbodies 


The CLOUDY setup for our models with blackbody SEDs are nearly identical to 
those with SSPs. To gain insight into the requirements for the two-photon contin- 
uum to dominate over the ionizing SED, we initialize hydrogen-only gas at constant 
temperature (measured from [O 11] 44363/5007) and systematically vary the den- 
sity and blackbody temperature (Figure 8). A weak UV turnover begins to appear 
at temperatures > 65,000 K, which is hotter than a typical O star. A strong UV 
turnover requires at least T > 90,000 K. The two-photon continuum is also sensitive 
to gas density because, at high densities, /-changing collisions will suppress the two- 
photon emission relative to Lya. This is seen in the bottom panel of Figure 8 where 
we vary gas density for a fixed blackbody temperature of 100,000 K. The UV turnover 
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Parameter Value 
n 103 em? 
Ton 105-5 K 
U 0.0 
Teas 104-3 K 
Zo/Zo —1.1 
Zc/Zo —1.6 
032 37.1 
He 11 leakage 0.22 


Table 3 Parameters for the 
optimized blackbody model. 


is strongly suppressed at n > 103 cm~’. We emphasize that the details of this cal- 
culation are sensitive to the chosen column density at which the model is truncated 
(which in our case is set by the velocity offset of Lya). Furthermore, there exist sig- 
nificant differences in atomic data predictions for the strength of l-changing collisions 
as a function of temperature [55]. 

To optimize the blackbody model fit to GS-NDG-9422, we begin with the param- 
eters of the best fit BPASS model and update ionization parameter, gas density, and 
blackbody temperature to simultaneously reproduce the emission line ratios and con- 
tinuum shape. We allow for both density and ionization bounded nebulae by adding 
an additional stopping criteria to reproduce the measured [O 11] \\3727/[O 111] 45007 
(O32) ratio. This stopping criteria often supersedes the neutral gas column stopping 
criteria used above. The O32 ratio and the gas temperature are allowed to vary within 
their observational uncertainties. As a final step, we add an additional parameter to 
extinguish the ionizing flux at E > 4 Rydberg to control the strength of the He 1 
emission, varying the leakage fraction of He l-ionizing photons from the blackbody 
between 0 and 1 (coloured lines in second panel of Figure 3). The parameters for our 
optimized model are reported in Table 3. We note that all our models predict strong 
Mg 11 A2796, 2803 emission at Acte © 1.95 um. Since this doublet is resonant and 
fese(Lya) < 100%, this feature is likely to be strongly suppressed. Alternatively, its 
absence in GS-NDG-9422 could simply indicate sub-solar Mg/O [56]. 


Models with Wolf-Rayet stars 


Adopting parameters from the optimised blackbody model but removing the He 1 
leakage parameter, we replaced the blackbody SED with a theoretical low-metallicity 
Wolf-Rayet star spectrum. We consider the WNL-H40 grid of the PoWR Wolf-Rayet 
models [24] with Z = 0.07 Zo, similar to that derived for GS-NDG-9422. We identify 
model 13-10 as most similar to our optimised blackbody, with T = 100,000 K and L = 
1053 Lo. No further optimisation is performed and the resulting spectrum is shown as 
the dashed orange line in Figure 3, nearly identical to the optimised blackbody model. 


Models with individual Population IIT stars 


We repeat the exercise described for our Wolf-Rayet models, this time replacing the 
blackbody SED with three Pop. III star models from [22] with effective temperatures 
between 95,000 K and 99,000 K. 
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Models with high-mass X-ray binaries stars 


Following [19, 57], we model the X-ray sources using a modified blackbody spectrum 
[58]. We assume black hole masses in the range 6 — 25 Mọ and disk radii between 
10° — 10 gravitational radii. In order for a UV turnover to appear, the ionizing photon 
output from the XRBs must dominate over the stellar population so that the nebular 
continuum can outshine the stars, implying a very high ratio of X-ray luminosity to star 
formation rate. In Figure 9, we show our XRB model for 25 Mọ black holes optimised 
to reproduce the continuum shape of GS-NDG-9422. This model significantly over- 
predicts the strength of the He 11 \1640 and \4686 lines. Varying the black hole masses 
does not resolve this issue. We therefore conclude that XRBs are not the dominant 
ionizing source in GS-NDG-9422. 


On the lack of Wolf-Rayet wind features 


Although Wolf-Rayet stars are typically associated with broad emission lines due to 
their high-velocity stellar winds, the absence of these features in GS-NDG-9422 could 
simply be the result of the extremely bright nebular component outshining these wind 
features, as predicted by our photoionization models. Furthermore, in the absence of 
iron, which is the dominant source of stellar atmospheric opacity, wind speeds can 
drop below 500 km s~!, even at solar oxygen abundance [59]. Hence, high-redshift 
Wolf-Rayet dominated galaxies may not show broad He 11 41640 and He 11 \4686 lines 
[60]. 


Ionizing flux required to power HG 


The hydrogen ionzing photon luminosity, Q, can be measured from the HO flux as: 


Aer D? Ing QB 
(1 ES fesc hung asi? 


Qe (1) 


Where D, = 58456.2 Mpc is the luminosity distance, Ing is the measured H£ flux, 
h is Planck’s constant, na is the frequency of the HO transition, ag is the total 
case B recombination rate, Lee is the escape fraction of ionzing photons, and oe is 
the effective HG recombination rate. To calculate a luminosity distance, we adopt a 
cosmology with Ho = 67.31 km s7! Mpc™! and Qm = 0.315 [61]. The escape fraction 
is derived from our best fitting blackbody photoionization model to be 7.3% and the 
recombination rates are evaluated at 20,000 K [62]. 


IMF excess 


To calculate the excess number of massive stars we re-run our best Wolf-Rayet model 
and the low-metallicity massive star model with Tog = 97,352 K while simultaneously 
adding a second component BPASS model. In the case of low-metallicity massive 
stars, we assume the BPASS model has an age of 3 Myr (the approximate lifetime of 
massive stars), while for the Wolf-Rayet models, we assume the progenitor stars are 
lower mass (50 Ma) and live for slightly longer (5 Myr). We assume the BPASS stellar 
population has the same metallicity as the gas. We then progressively increase the 
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ionization parameter of the second population until the UV turnover from the two- 
photon continuum begins to weaken. We calculate the stellar mass of a BPASS SSP 
that can be present per hot star as 


Ussp Qstar 
Ustar Qssp 


Mssp = (2) 


where Qstar = 10493 y s71 for Wolf-Rayet stars or 1049-98 y s~! for low-metallicity 
massive stars with Teg = 97,352 K. Qssp = 1046:73 or 1046-03 y s-1Mzt at 3 Myr and 
5 Myr, respectively. In both cases pase = = 6.3%. Therefore, Mssp = 112 or 137 Mo for 
low-metallicity massive stars and Wolf-Rayet stars, respectively. Based on a BPASS 
SSP of 112 Mo, one would expect on average 0.029 stars of > 100 Mo while for an 
SSP of 137 Mo, the expectation is 0.105 stars of > 50 Mo. In other words, there is a 
35x excess in the number of massive stars in the case of our low-metallicity massive 
star model, or a 9.5x excess for our Wolf-Rayet model. 

The results of this calculation are sensitive to the chosen underlying SSP, the pro- 
genitor masses of Wolf-Rayet stars, and the mass of the low-metallicity massive stars. 
For SSPs that have intrinsically lower Zon, the excess drops as Qssp appears in the 
denominator of Figure 2. If the ages at which the exotic stars evolve off the main 
sequence is longer than what we have assumed here, the excess will also decrease. 
Nevertheless, we highlight that the excess number of massive stars we have calcu- 
lated represents a conservative lower limit because we have maximized the allowable 
contribution from the SSP in the model. 

These values allow us to calculate the maximum mass of the star clusters that 
formed in GS-NDG-9422. Calculating the number of stars based on the HO luminosity, 
GS-NDG-9422 would be powered by ~ 17,000 very metal-poor stars of ~ 100 Mo with 
Teg = 97,000 K, or ~ 70,000 low-metallicity Wolf-Rayet stars. Based on the above 
calculations, this implies a maximum star cluster mass of 10°°° Mo or 1071? Mo for 
the low-metallicity massive star or Wolf-Rayet star cases, respectively. 
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Fig. 8 (Top) Normalized spectrum of hydrogen-only gas with n = 103 cm~? irradiated by blackbod- 
ies of different temperatures (as given in the colour bar). The magenta line represents a blackbody 
temperature of 100,000 K. (Bottom) Normalized spectrum of hydrogen-only gas irradiated by a 
blackbody with T = 10° K at different gas densities (as given in the colour bar). The magenta line 
represents a density of 103 em". 
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Fig. 9 JWST NIRSpec PRISM spectrum of GS-NDG-9422 (black) compared with photoionization 
models that include high-mass X-ray binaries with a black hole mass of 25 Mọ (magenta). The solid 
and dashed magenta lines indicate models with different accretion disk radii. 
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